function [ subdiv pick123 pick4 ] = loop_subdivmatrix(valency)

% LOOP_SUBDIVMATRIX  Returns the Loop subdivision matrix
%   and 'picking matrices' for the given valency

% --
% Memoize for performance. All valencies must lie between 3 and max_val
% --
max_val = 20;
persistent all_subdiv ord_pick123 ord_pick4;
persistent ext_pick123_val3 ext_pick123_stub;

if isempty(all_subdiv)
    all_subdiv       = cell(max_val - 2);
end
if isempty(ord_pick123)
    ord_pick123      = zeros(12, 18, 3);
end
if isempty(ord_pick4)
    ord_pick4        = zeros(12, 18);
end
if isempty(ext_pick123_val3)
    ext_pick123_val3 = zeros(12, 15, 3);
end
if isempty(ext_pick123_stub)
    ext_pick123_stub = zeros(12, 16, 3);
end

if valency == 6 && ~isempty(all_subdiv{4})
    subdiv  = all_subdiv{4};
    pick123 = ord_pick123;
    pick4   = ord_pick4;
    return;
elseif ~isempty(all_subdiv{valency - 2})
    subdiv  = ext_subdiv{valency - 2};
    if valency == 3
        pick123 = ext_pick123_val3;
    else
        pick123 = [ ext_pick123_stub zeros(12, valency - 4, 3) ];
    end
    pick4 = [ zeros(6 + valency, 6) eye(6 + valency) ];
    return;
end

% --
% Normal function
% --

pick123 = zeros(12, 12 + valency, 3);

if valency == 6
    subdiv = [ 0.375  0.125  0      0.125  0.375  0      0      0      0      0      0      0      ;
               0.125  0.375  0      0      0.375  0.125  0      0      0      0      0      0      ;
               0      0.375  0.125  0      0.125  0.375  0      0      0      0      0      0      ;
               0      0.125  0.375  0      0      0.375  0.125  0      0      0      0      0      ;
               0.125  0      0      0.375  0.375  0      0      0.125  0      0      0      0      ;
               0.0625 0.0625 0      0.0625 0.625  0.0625 0      0.0625 0.0625 0      0      0      ;
               0      0.125  0      0      0.375  0.375  0      0      0.125  0      0      0      ;
               0      0.0625 0.0625 0      0.0625 0.625  0.0625 0      0.0625 0.0625 0      0      ;
               0      0      0.125  0      0      0.375  0.375  0      0      0.125  0      0      ;
               0      0      0      0.125  0.375  0      0      0.375  0.125  0      0      0      ;
               0      0      0      0      0.375  0.125  0      0.125  0.375  0      0      0      ;
               0      0      0      0      0.125  0.375  0      0      0.375  0.125  0      0      ;
               0      0      0      0      0      0.375  0.125  0      0.125  0.375  0      0      ;
               0      0      0      0      0.125  0      0      0.375  0.375  0      0.125  0      ;
               0      0      0      0      0.0625 0.0625 0      0.0625 0.625  0.0625 0.0625 0.0625 ;
               0      0      0      0      0      0.125  0      0      0.375  0.375  0      0.125  ;
               0      0      0      0      0      0      0      0.125  0.375  0      0.375  0.125  ;
               0      0      0      0      0      0      0      0      0.375  0.125  0.125  0.375  ];

    pick123(:, :, 1) = [ 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ;
                         0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ;
                         0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ;
                         0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 ;
                         0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 ;
                         0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 ;
                         0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 ;
                         0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 ;
                         0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 ;
                         0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 ;
                         0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 ;
                         0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 ];
                     
    pick123(:, :, 2) = [ 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 ;
                         0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 ;
                         0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 ;
                         0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 ;
                         0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 ;
                         0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 ;
                         0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 ;
                         0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 ;
                         0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 ;
                         0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 ;
                         0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ;
                         0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ];
                     
    pick123(:, :, 3) = [ 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ;
                         0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ;
                         0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ;
                         0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 ;
                         0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 ;
                         0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 ;
                         0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 ;
                         0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 ;
                         0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 ;
                         0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 ;
                         0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 ;
                         0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 ];
                     
    pick4            = [ 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 ;
                         0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 ;
                         0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 ;
                         0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 ;
                         0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 ;
                         0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 ;
                         0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 ;
                         0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 ;
                         0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 ;
                         0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 ;
                         0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 ;
                         0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 ];
                     
    all_subdiv{valency - 2} = subdiv;
    ord_pick123 = pick123;
    ord_pick4   = pick4;
else    
    alpha = (0.375 + cos(2 * pi / valency) / 4) ^ 2 + 0.375;
    delta = (1 - alpha) / valency;
    subdiv = zeros(6 + valency);

    subdiv(1, 1:9) = [ 1 1 0 1 10  1 0 1 1 ];
    subdiv(2, 1:9) = [ 0 2 0 0  6  6 0 0 2 ];
    subdiv(3, 1:9) = [ 0 1 1 0  1 10 1 0 1 ];
    if valency == 3, subdiv(3, 8) = 1; else subdiv(3, 10) = 1; end
    subdiv(4, 1:9) = [ 0 0 0 2  6  0 0 6 2 ];
    subdiv(5, 1:9) = [ 0 0 0 0  6  2 0 2 6 ];
    subdiv(6, 1:9) = [ 0 0 0 0  2  6 0 0 6 ];
    if valency == 3, subdiv(6, 8) = 2; else subdiv(6, 10) = 2; end
    subdiv(7, 1:9) = [ 0 0 0 0  0  6 2 0 2 ];
    if valency == 3, subdiv(7, 8) = 6; else subdiv(7, 10) = 6; end
    subdiv(8, 1:9) = [ 0 0 0 0  2  0 0 6 6 ];
    if valency > 4
        subdiv(8, 11) = 2;
    elseif valency == 4
        subdiv(8, 10) = 2;
    else
        assert(valency == 3);
        subdiv(8, 6) = 2;
    end
    subdiv(9, 1:9)  = [ 0 0 0 0 delta delta 0 delta alpha ];
    subdiv(9, 10:end) = delta;

    if valency > 3
        subdiv(10, 6:10) = [ 2 0 0 6 6 ];

        if valency == 4
            subdiv(10, 8) = 2;
        else
            subdiv(10, end) = 2;
        end
    end

    if valency > 4
        subdiv(11, 8:11) = [ 2 6 0 6 ];

        if valency == 5
            subdiv(11, 10) = 2;
        else
            subdiv(11, 12) = 2;
        end
    end

    if valency > 5
        for pt = 12:valency + 6
            subdiv(pt, 9) = 6;
            subdiv(pt, pt - 1:pt) = [ 2 6 ];

            if pt == valency + 6
                subdiv(pt, 10) = 2;
            else
                subdiv(pt, pt + 1) = 2;
            end
        end
    end
    
    extras = zeros(6, 6 + valency);

    extras(:, 1:8) = [ 2 0 0 6 6 0 0 2 ;
                       6 2 0 2 6 0 0 0 ;
                       2 6 0 0 6 2 0 0 ;
                       0 6 2 0 2 6 0 0 ;
                       0 2 6 0 0 6 2 0 ;
                       0 0 2 0 0 6 6 0 ];

    if valency == 3
        extras(6, 8) = 2;
    else
        extras(6, 10) = 2;
    end
    
    extras = extras ./ 16;
    
    subdiv = subdiv ./ repmat(sum(subdiv, 2), 1, 6 + valency);
    subdiv = [ extras ; subdiv ];
    
    ext_subdiv{valency - 2} = subdiv;
    
    pick4 = [ zeros(6 + valency, 6) eye(6 + valency) ];
    
    pick123(:, 1:15, 1) = [ 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 ;
                            0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 ;
                            0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 ;
                            1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ;
                            0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 ;
                            0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 ;
                            0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 ;
                            0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 ;
                            0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 ;
                            0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 ;
                            0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 ;
                            0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 ];
                        
    if valency == 3
        pick123(:, 1:15, 2) = [ 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 ;
                                0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 ;
                                0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 ;
                                0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 ;
                                0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 ;
                                0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 ;
                                0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 ;
                                0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 ;
                                0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 ;
                                0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 ;
                                0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 ;
                                0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 ];
                            
        pick123(:, 1:15, 3) = [ 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 ;
                                0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 ;
                                0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 ;
                                0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 ;
                                0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 ;
                                0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 ;
                                0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 ;
                                0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 ;
                                0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 ;
                                0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 ;
                                0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 ;
                                0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 ];

        ext_pick123_val3 = pick123;
    else
        pick123(:, 1:16, 2) = [ 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 ;
                                0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 ;
                                0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 ;
                                0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 ;
                                0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 ;
                                0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 ;
                                0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 ;
                                0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 ;
                                0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 ;
                                0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 ;
                                0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 ;
                                0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 ];
                            
        pick123(:, 1:16, 3) = [ 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 ;
                                0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 ;
                                0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 ;
                                0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 ;
                                0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 ;
                                0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 ;
                                0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 ;
                                0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 ;
                                0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 ;
                                0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 ;
                                0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 ;
                                0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 ];
                            
        ext_pick123_stub = pick123(:, 1:16, :);
    end
end

end